//*********************************************
//  This is Geant4 Template
//                                  author:Qian
//

#include "G4AnalysisManager.hh"
#include "MyAnalysisManager.hh"
#include "Verbose.hh"

#include "G4AntiNeutrinoE.hh"
#include "G4AntiNeutrinoMu.hh"
#include "G4AntiNeutrinoTau.hh"

#include "G4ProcessType.hh"
//fNotDefined, fTransportation, fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
//fDecay, fGeneral, fParameterisation, fUserDefined, fParallel, fPhonon, fUCN
#include "G4DecayProcessType.hh"
//DECAY, DECAY_WithSpin, DECAY_PionMakeSpin, DECAY_Radioactive, DECAY_Unknown, DECAY_External
#include "G4HadronicProcessType.hh"
//fHadronElastic, fHadronInelastic, fCapture, fFission, fHadronAtRest, fLeptonAtRest, fChargeExchange, fRadioactiveDecay
#include "G4TransportationProcessType.hh"
//TRANSPORTATION, COUPLED_TRANSPORTATION, STEP_LIMITER, USER_SPECIAL_CUTS, NEUTRON_KILLER
#include "G4StepStatus.hh"
//fWorldBoundary, fGeomBoundary, fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc, fUserDefinedLimit, fExclusivelyForcedProc, fUndefined
#include "G4TrackStatus.hh"
//fAlive, fStopButAlive, fStopAndKill, fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent
#include "G4Track.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4RegionStore.hh"
#include "G4Region.hh"

#include "MyDetectorConstruction.hh"
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wshadow"
#include "TH1F.h"
#pragma GCC diagnostic pop

#include "MySingleCsIAnalysis.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

MySingleCsIAnalysis::MySingleCsIAnalysis()
{
    SetFileName("g4output.root");

    //-------
    //#ANALYSIS 1. 初始化变量

    fTrkLen = 0;
    fEdeps.clear();
    fHitsX.clear();
    fHitsY.clear();
    fHitsZ.clear();
    Eng_Ei.clear();
    Eng_MSC.clear();
    Eng_SC.clear();
    Eng_Pho.clear();
    num=0;
    Nhits=0;
    Nhits1=0;
    Nhits2=0;
    
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

MySingleCsIAnalysis::~MySingleCsIAnalysis()
{
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void MySingleCsIAnalysis::BeginOfRunAction()
{
    
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MySingleCsIAnalysis::BeginOfRunAction()" << G4endl;

    auto analysisManager = G4AnalysisManager::Instance();

    // Default settings
    analysisManager->SetNtupleMerging(true);
    // Note: merging ntuples is available only with Root output

    analysisManager->SetVerboseLevel(1);
    analysisManager->OpenFile(fFileName);

    //-------
    //#ANALYSIS 2. 定义Ntuple结构

    // Creating 1D histograms
    //analysisManager->SetFirstHistoId(1);
    //analysisManager->CreateH1("phEng", "photon energy", 50, 0., 100); // h1 Id = 0

    // Creating 2D histograms
    //analysisManager->CreateH2("HitOnAnode", "Cherenkov photon hits on the anodes", // h2 Id = 0
    //                          50, -1000., 1000, 50, -1000., 1000.);

    // Creating ntuple
    //
    analysisManager->SetFirstNtupleId(1);

    analysisManager->CreateNtuple("Hits", "Hits"); // ntuple Id = 1
    analysisManager->CreateNtupleDColumn("Nhits");//Nhits1+Nhits2 every event
    analysisManager->CreateNtupleDColumn("Nhits1");
    analysisManager->CreateNtupleDColumn("Nhits2");
    analysisManager->CreateNtupleDColumn("num");//the number of event
    analysisManager->CreateNtupleDColumn("Sum");
    analysisManager->CreateNtupleDColumn("e_eng");
    analysisManager->CreateNtupleDColumn("opt_eng");
    analysisManager->CreateNtupleDColumn("gamma_dep");
    // analysisManager->CreateNtupleDColumn("Gamma_Euger");
    // analysisManager->CreateNtupleDColumn("Gamma_Fluo");
    // analysisManager->CreateNtupleDColumn("Gamma_PIXE");
    // analysisManager->CreateNtupleDColumn("E_PIXE");
    // analysisManager->CreateNtupleDColumn("E_Fluo");
    // analysisManager->CreateNtupleDColumn("E_Euger");
    analysisManager->CreateNtupleDColumn("Eng_Ele_lose");
    analysisManager->CreateNtupleDColumn("num_mod_0");
    analysisManager->CreateNtupleDColumn("num_mod_1");
    analysisManager->CreateNtupleDColumn("mod");
    analysisManager->CreateNtupleDColumn("Eng_e_sum");
    analysisManager->CreateNtupleDColumn("Eng_g_sum");
    analysisManager->CreateNtupleDColumn("Eng_g_0_dop_sum");
    // analysisManager->CreateNtupleDColumn("Eng_Pho_all");
    analysisManager->CreateNtupleDColumn("MSC");
    analysisManager->CreateNtupleDColumn("EIoni");
    analysisManager->CreateNtupleDColumn("SCintillation");
    analysisManager->CreateNtupleDColumn("Energy");
    analysisManager->CreateNtupleDColumn("E_num");
    analysisManager->CreateNtupleDColumn("Eng_e_1",Eng_e_1);
    analysisManager->CreateNtupleDColumn("Eng_g_1",Eng_g_1);
    analysisManager->CreateNtupleDColumn("Eng_e_1_dop",Eng_e_1_dop);
    analysisManager->CreateNtupleDColumn("Eng_g_1_dop",Eng_g_1_dop);
    analysisManager->CreateNtupleDColumn("Eng_g_0_dop",Eng_g_0_dop);
    analysisManager->CreateNtupleDColumn("Eng_MSC",Eng_MSC);//deposit energy
    analysisManager->CreateNtupleDColumn("Eng_Ei",Eng_Ei);
    analysisManager->CreateNtupleDColumn("Eng_SC",Eng_SC);
    analysisManager->CreateNtupleDColumn("Eng_Pho",Eng_Pho);
    analysisManager->FinishNtuple();

    // analysisManager->CreateNtuple("Photon", "Photon"); // ntuple Id = 1
    // analysisManager->CreateNtupleDColumn("Eng_Pho",Eng_Pho);
    // analysisManager->FinishNtuple();

    return;
}

void MySingleCsIAnalysis::EndOfRunAction()
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MySingleCsIAnalysis::EndOfRunAction()" << G4endl;

    //-------
    //#ANALYSIS 6. 在Run结束的时候将ntuple保存到文件
    
    auto analysisManager = G4AnalysisManager::Instance();

    analysisManager->Write();
    analysisManager->CloseFile();

    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void MySingleCsIAnalysis::BeginOfEventAction(const G4Event *)
{
    if (!active)
        return;

    if (verbose > 1)
        G4cout << "====>MySingleCsIAnalysis::BeginOfEventAction()" << G4endl;

    //-------
    //#ANALYSIS 3. 初始化Event开始的参数
    fTrkLen = 0;
    fEdeps.clear();
    fHitsX.clear();
    fHitsY.clear();
    fHitsZ.clear();
    Eng_Ei.clear();
    Eng_SC.clear();
    Eng_MSC.clear();
    Eng_Pho.clear();
    Eng_e_1_dop.clear();
    Eng_g_1_dop.clear();
    Eng_e_1.clear();
    Eng_g_1.clear();
    Eng_g_0_dop.clear();
    num_mod_0=0;
    num_mod_1=0;
    mod=0;
    Eng_e_sum=0;
    Eng_g_sum=0;
    Eng_g_0_dop_sum=0;

    Eng_Pho_all=0;
    Nhits=0;
    Nhits1=0;
    Nhits2=0;
    MSC=0;
    EIoni=0;
    SCintillation=0;
    Sum=0;
    num=0;
    Energy=0;
    E_num=0;
    Select=0;
    Eng_Ele_Lose=0;

    //model:
    Gamma_Euger=0;
    Gamma_Fluo=0;
    Gamma_PIXE=0;
    E_PIXE=0;
    E_Fluo=0;
    E_Euger=0;

    //energy
    e_eng=0;
    gamma_dep=0;
    opt_eng=0;
    return;
}

void MySingleCsIAnalysis::EndOfEventAction(const G4Event *)
{
    if (!active)
        return;

    if (verbose > 1)
        G4cout << "====>MySingleCsIAnalysis::EndOfEventAction()" << G4endl;

    //-------
    //#ANALYSIS 5. 在Event结束的时候将数据保存到ntuple
    Nhits=Nhits1+Nhits2;
    //num=num+1;
    // G4cout<<"Nhits"<<Nhits<<G4endl;
    // G4cout<<"num"<<num<<G4endl;
    // G4cout<<"energy:"<<Energy<<G4endl;
    auto analysisManager = G4AnalysisManager::Instance();
    analysisManager->FillNtupleDColumn(1, 0, Nhits);
    analysisManager->FillNtupleDColumn(1, 1, Nhits1);
    analysisManager->FillNtupleDColumn(1, 2, Nhits2);
    analysisManager->FillNtupleDColumn(1, 3, num);
    analysisManager->FillNtupleDColumn(1, 4, Sum);
    analysisManager->FillNtupleDColumn(1, 5, e_eng);
    analysisManager->FillNtupleDColumn(1, 6, opt_eng);
    analysisManager->FillNtupleDColumn(1, 7, gamma_dep);
    analysisManager->FillNtupleDColumn(1, 8, Eng_Ele_Lose);
    analysisManager->FillNtupleDColumn(1, 9, num_mod_0);
    analysisManager->FillNtupleDColumn(1, 10, num_mod_1);
    analysisManager->FillNtupleDColumn(1, 11, mod);
    analysisManager->FillNtupleDColumn(1, 12, Eng_e_sum);
    analysisManager->FillNtupleDColumn(1, 13, Eng_g_sum);
    analysisManager->FillNtupleDColumn(1, 14, Eng_g_0_dop_sum);
    // analysisManager->FillNtupleDColumn(1, 6, Eng_Pho_all);
    // analysisManager->FillNtupleDColumn(1, 4, MSC);
    // analysisManager->FillNtupleDColumn(1, 5, EIoni);
    // analysisManager->FillNtupleDColumn(1, 6, SCintillation);
    
    // analysisManager->FillNtupleDColumn(1, 8, Energy);
    // analysisManager->FillNtupleDColumn(1, 9, E_num);
    // analysisManager->FillNtupleDColumn(1, 10, Eng_Pho_all);
    //analysisManager->FillNtupleDColumn(1, 4, Eng);
    analysisManager->AddNtupleRow(1);

    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ClassificationOfNewTrack MySingleCsIAnalysis::ClassifyNewTrack(const G4Track *aTrack)
{
    if (!active)
        return fUrgent;

    if (verbose > 2)
        G4cout << "====>MySingleCsIAnalysis::ClassifyNewTrack()" << G4endl;
    if (aTrack->GetParentID() == 0) { return fUrgent; }

    if(!fIDdefined) {
    fIDdefined = true;
    fPhotoGamma = G4PhysicsModelCatalog::GetModelIndex("phot_fluo");
    fComptGamma = G4PhysicsModelCatalog::GetModelIndex("compt_fluo");
    fPhotoAuger = G4PhysicsModelCatalog::GetModelIndex("phot_auger");
    fComptAuger = G4PhysicsModelCatalog::GetModelIndex("compt_auger");
    fPixeGamma = G4PhysicsModelCatalog::GetModelIndex("gammaPIXE");
    fPixeAuger = G4PhysicsModelCatalog::GetModelIndex("e-PIXE");
    fElectronDNAGamma =
     G4PhysicsModelCatalog::GetModelIndex("e-_G4DNAIonisation_fluo");
    fElectronDNAAuger =
     G4PhysicsModelCatalog::GetModelIndex("e-_G4DNAIonisation_auger");
    fProtonDNAGamma =
     G4PhysicsModelCatalog::GetModelIndex("proton_G4DNAIonisation_fluo");
    fProtonDNAAuger =
     G4PhysicsModelCatalog::GetModelIndex("proton_G4DNAIonisation_auger");
    fHydrogenDNAGamma =
     G4PhysicsModelCatalog::GetModelIndex("hydrogen_G4DNAIonisation_fluo");
    fHydrogenDNAAuger =
     G4PhysicsModelCatalog::GetModelIndex("hydrogen_G4DNAIonisation_auger");
    fAlphaDNAGamma =
     G4PhysicsModelCatalog::GetModelIndex("alpha_G4DNAIonisation_fluo");
    fAlphaDNAAuger =
     G4PhysicsModelCatalog::GetModelIndex("alpha_G4DNAIonisation_auger");
    fAlphaPlusDNAGamma =
     G4PhysicsModelCatalog::GetModelIndex("alpha+_G4DNAIonisation_fluo");
    fAlphaPlusDNAAuger =
     G4PhysicsModelCatalog::GetModelIndex("alpha+_G4DNAIonisation_auger");
    fHeliumDNAGamma =
     G4PhysicsModelCatalog::GetModelIndex("helium_G4DNAIonisation_fluo");
    fHeliumDNAAuger =
     G4PhysicsModelCatalog::GetModelIndex("helium_G4DNAIonisation_auger");
    fGenericIonDNAGamma =
     G4PhysicsModelCatalog::GetModelIndex("GenericIon_G4DNAIonisation_fluo");    
     fGenericIonDNAAuger =
     G4PhysicsModelCatalog::GetModelIndex("GenericIon_G4DNAIonisation_auger");
  }
    
    //-------
    //#ANALYSIS 4.1 在生成新Track的时候保存相应数据
        G4ParticleDefinition *particle = aTrack->GetDefinition();
        G4String name_track = particle->GetParticleName();
        G4double Eng=aTrack->GetTotalEnergy();
        
        if(name_track=="opticalphoton"){
            num=num+1;
            opt_eng=Eng+opt_eng;
        }
        
        
        
        
        // if(name_track=="e-"){
        //     e_eng=Eng+e_eng;
        // }
            
        
        //G4cout<<"model: "<<idx<<G4endl;
        
        // const G4ParticleDefinition *part = aTrack->GetDefinition();
        // G4String partType = part->GetParticleName();
        // const G4String procName = aTrack->GetCreatorProcess()->GetProcessName();
        // G4cout << "粒子类型：" << partType << " – "<< "过程类型：" << procName << "粒子能量：" << aTrack->GetKineticEnergy() << G4endl;
        //G4cout<<"Model"<<creat_pro<<G4endl;
        //G4cout<<"phot_fluo:"<<G4PhysicsModelCatalog::GetIndex("phot_fluo")<<G4endl;
        //G4cout<<"e fluo "<<G4PhysicsModelCatalog::GetIndex("e-_G4DNAIonisation_fluo")<<G4endl;
        //开始判断：
        // if(idx==fPhotoGamma)
        // if(idx==fPhotoAuger)
        //Eng_Pho.push_back(aTrack->GetTotalEnergy());
        if(Select==1){
            //if(idx==fPhotoGamma)G4cout<<"this is fluo"<<G4endl;
            
            
            //Eng_Pho_all=Ene_Pho+Eng_Pho_all;
            
        }
    
    //if (aTrack->GetParticleDefinition() == G4OpticalPhoton::G4OpticalPhoton())
    //{
    //    return fKill;
    //
    //auto analysisManager = G4AnalysisManager::Instance();
   // analysisManager->FillNtupleDColumn(4, 0, aTrack->GetKineticEnergy());
    //analysisManager->AddNtupleRow(4);

    return fUrgent;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MySingleCsIAnalysis::PreTrackingAction(const G4Track *aTrack)
{
    if (!active)
        return;

    if (verbose > 2)
        G4cout << "====>MySingleCsIAnalysis::PreTrackingAction()" << G4endl;

        //G4double ene_pho = aTrack->GetKineticEnergy();
        
        //G4ThreeVector trkPos = aTrack->GetPosition();
        //G4cout<<"the trkPos is"<<trkPos<<G4endl;
        //G4TrackStatus trkStatus = aTrack->GetTrackStatus(); //检查track的状态，常用fStopAndKill判断
        //G4cout<<"the trkStatus is "<<trkStatus<<G4endl;  
       
            //还可根据坐标判断是否在Region内
        /*const G4ThreeVector& position = aTrack->GetPosition();

        G4Region *region = G4RegionStore::GetInstance()->GetRegion("SiPMVol");
        int N =  region->GetNumberOfRootVolumes();
        std::vector<G4LogicalVolume*>::iterator it_logicalVolumeInRegion = region->GetRootLogicalVolumeIterator();

        bool inside_target = false;

        for(int i = 0; i < N ; i++, it_logicalVolumeInRegion++)
        {
            EInside test_status = (*it_logicalVolumeInRegion)->GetSolid()->Inside(position) ;
            if(test_status == kInside)
            {
                //...
                break;
            }
            else if (test_status == kSurface)
            {
                //...
                break;
            }
        }*/

        //...
        // G4ParticleDefinition *particle = aTrack->GetDefinition();
        // G4String name_track = particle->GetParticleName();
        // if(Select==1){
        //     if(name_track=="opticalphoton")
        //     Eng_Pho.push_back(aTrack->GetTotalEnergy());
        // }
    //通常用法为：
    //  保存某类粒子的相关信息(更多的时候是母粒子的次级粒子)，如产生点位置，产生点时间，物理过程，粒子种类等。
    //  特别是衰变的粒子或打靶实验
    /*
    {
        G4ParticleDefinition *particle = aTrack->GetDefinition();
        G4String name = particle->GetParticleName();
        G4int pid = particle->GetPDGEncoding();
        G4int Z = particle->GetAtomicNumber();
        G4int A = particle->GetAtomicMass();
        G4double charge = particle->GetPDGCharge();

        G4int ID = aTrack->GetTrackID();
        G4double energy = aTrack->GetKineticEnergy();
        
        G4double t0 = aTrack->GetGlobalTime();
        G4double weight = aTrack->GetWeight();
        G4int processType = aTrack->GetCreatorProcess()->GetProcessSubType(); //查看ProcessType的头文件定义
        G4String processName = aTrack->GetCreatorProcess()->GetProcessName();

        //还可根据坐标判断是否在Region内
        const G4ThreeVector& position = atrack->GetPosition();

        G4Region *region = G4RegionStore::GetInstance()->GetRegion("regName");
        int N =  region->GetNumberOfRootVolumes();
        std::vector<G4LogicalVolume*>::iterator it_logicalVolumeInRegion = region->GetRootLogicalVolumeIterator();

        bool inside_target = false;

        for(int i = 0; i < N ; i++, it_logicalVolumeInRegion++)
        {
            EInside test_status = (*it_logicalVolumeInRegion)->GetSolid()->Inside(position) ;
            if(test_status == kInside)
            {
                //...
                break;
            }
            else if (test_status == kSurface)
            {
                //...
                break;
            }
        }

        //...
    }
    */

    //-------
    //#ANALYSIS 4.2 在Tracking产生的时候保存相应数据
    G4ParticleDefinition *particle = aTrack->GetDefinition();
    G4String name_track = particle->GetParticleName();
    G4int parentID = aTrack->GetParentID();
    G4int idx = aTrack->GetCreatorModelID();
    G4String creat_pro=aTrack->GetCreatorModelName();
    G4double Kine=aTrack->GetKineticEnergy();
    G4String processName;
    //if(parentID==1)processName = aTrack->GetCreatorProcess()->GetProcessName();
    //G4String processName = aTrack->GetCreatorProcess()->GetProcessName();
    if(parentID!=1&&parentID!=0)processName="sha";
    if(idx==32)mod++;
        if (parentID==1&&name_track=="e-")
        {
            processName = aTrack->GetCreatorProcess()->GetProcessName();
            if(processName=="pol-phot"){
            num_mod_0++;
            Eng_e_1.push_back(aTrack->GetKineticEnergy());
            Eng_e_sum=Eng_e_sum+Kine;
            }
        }
        if(parentID==1&&name_track=="gamma"){
            
            if(processName=="pol-phot"){
            num_mod_0++;
            Eng_g_1.push_back(aTrack->GetKineticEnergy());
            Eng_g_sum=Eng_g_sum+Kine;}
        }
    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MySingleCsIAnalysis::PostTrackingAction(const G4Track *)
{
    if (!active)
        return;

    if (verbose > 2)
        G4cout << "====>MySingleCsIAnalysis::PostTrackingAction()" << G4endl;
        //G4ThreeVector trkPos = aTrack->GetPosition();
        //G4cout<<"the trkPos is"<<trkPos<<G4endl;
        //G4TrackStatus trkStatus = aTrack->GetTrackStatus(); //检查track的状态，常用fStopAndKill判断
        //G4cout<<"the trkStatus is "<<trkStatus<<G4endl;
    //通常用法为：
    //  1. 根据粒子停止位置，计算能量泄漏等。
    //  2. 根据粒子停止时间，计算飞行时间等。
    //  3. 径迹总长度等信息
    /*
    {
        G4ParticleDefinition *particle = aTrack->GetDefinition();
        G4String name = particle->GetParticleName();
        G4int pid = particle->GetPDGEncoding();
        G4int Z = particle->GetAtomicNumber();
        G4int A = particle->GetAtomicMass();
        G4double charge = particle->GetPDGCharge();

        G4ThreeVector trkPos = aTrack->GetPosition();
        G4int ID = aTrack->GetTrackID();
        G4double energy = aTrack->GetKineticEnergy();
        G4double t1 = aTrack->GetGlobalTime();
        G4double weight = aTrack->GetWeight();
        G4int processType = aTrack->GetCreatorProcess()->GetProcessSubType(); //查看ProcessType的头文件定义
        G4String processName = aTrack->GetCreatorProcess()->GetProcessName();

        G4StepStatus stepStatus = aTrack->GetStep()->GetPostStepPoint()->GetStepStatus(); //检查最后一个step的状态, 常用fWorldBoundary判断leakage
        G4int nbSteps = aTrack->GetCurrentStepNumber();
        G4double Trleng = aTrack->GetTrackLength();

        G4TrackStatus trkStatus = aTrack->GetTrackStatus(); //检查track的状态，常用fStopAndKill判断
    }
    */

    //-------
    //#ANALYSIS 4.3 在Tracking终止的时候保存相应数据

    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MySingleCsIAnalysis::SteppingAction(const G4Step *aStep)
{
    if (!active)
        return;

    if (verbose > 2)
        G4cout << "====>MySingleCsIAnalysis::SteppingAction()" << G4endl;

    //---
    //1. 相关参数的获取

    //1.1 Track的相关参数
    const G4Track *aTrack = aStep->GetTrack();
    G4int parentID = aTrack->GetParentID();
    const G4ParticleDefinition *particle = aTrack->GetParticleDefinition();
    //G4int pdgID = particle->GetPDGEncoding();
    //G4int charge = particle->GetPDGCharge();
  
    //以下是G4Track常见的一些参数获取方法
    {
        /*
        G4int iStep = aTrack->GetCurrentStepNumber();
        G4int trackID = aTrack->GetTrackID();
        G4int parentID = aTrack->GetParentID();
        G4ParticleDefinition particle = aTrack->GetParticleDefinition();
        G4int pdgID = particle->GetPDGEncoding();
        G4String pdgName = particle->GetParticleName();
        G4cout << "-->" << iStep << ": tid=" << trackID << " " << pdgName << "(" << pdgID << ") parentID=" << parentID << G4endl;

        G4double energy = aTrack->GetKineticEnergy(); // kinetic energy at current step
        G4double engtot = aTrack->GetTotalEnergy();   // total energy (including m0)
        G4double velocity = aTrack->GetVelocity();    // mm/s
        G4ThreeVector moment = aTrack->GetMomentum();
        G4ThreeVector momDir = aTrack->GetMomentumDirection(); //unit vector
        G4ThreeVector polar = aTrack->GetPolarization();
        G4cout << "-->" << energy << " " << engtot << " " << velocity << " " << moment << " " << momDir << G4endl;

        //two different ways to get the process name
        G4String crtName = (parentID != 0) ? aTrack->GetCreatorModelName() : "none";
        G4String proName = (parentID != 0) ? aTrack->GetCreatorProcess()->GetProcessName() : "none";
        G4cout << "-->" << crtName << " " << proName << G4endl;

        G4ThreeVector trkPos = aTrack->GetPosition(); //current position of this track at the postStepPoint of this step
        G4double trkLen = aTrack->GetTrackLength(); //current track-length of this track at this postStepPoint of this step
        G4double gtime = aTrack->GetGlobalTime();
        G4double ltime = aTrack->GetLocalTime();
        G4double ptime = aTrack->GetProperTime();
        G4cout << "-->trkPos=" << trkPos << " " << trkLen << " " << gtime << " " << ltime << " " << ptime << G4endl;
        */
    }

    //1.2 Step的相关参数
    G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
    G4StepPoint *postStepPoint = aStep->GetPostStepPoint();
    auto *postVolume = postStepPoint->GetTouchableHandle()->GetVolume();
    auto *preVolume = preStepPoint->GetTouchableHandle()->GetVolume();
    if (postVolume == NULL)
        return;
    G4String PostVolume = postVolume->GetName();
    G4String PreVolume = preVolume->GetName();
    G4String name = particle->GetParticleName();
    
    //G4double preEngTot = preStepPoint->GetTotalEnergy();
    G4ThreeVector postPos = postStepPoint->GetPosition();
    //double z=postPos.z();
    // G4StepStatus states=postStepPoint->GetStepStatus();
    
    // auto *pVolume = postStepPoint->GetTouchableHandle()->GetVolume();
    // G4LogicalVolume *presentVolume = pVolume->GetLogicalVolume();
    if(name=="opticalphoton"){
        if(PostVolume=="PMTPhys1"&&PreVolume=="CsIPhys")
       // if(z<-140&&PostVolume=="PMTPhys1"&& aStep->IsFirstStepInVolume())
        {
            Nhits1=Nhits1+1;
            //G4cout<<"the first Nhits "<<Nhits1<<G4endl;
        }
        if(PostVolume=="BoxPhys"&&PreVolume=="CsIPhys")Nhits=Nhits+1;
        if(PostVolume=="CsIPhys"&&PreVolume=="PMTPhys1")Nhits1=Nhits1-1;

        //if(z>140&&PostVolume=="PMTPhys2"&& aStep->IsFirstStepInVolume())
        if(PostVolume=="PMTPhys2"&&PreVolume=="CsIPhys")
        {
            Nhits2=Nhits2+1;
            //G4cout<<"the second Nhits "<<Nhits2<<G4endl;
        }
        if(PostVolume=="CsIPhys"&&PreVolume=="PMTPhys2")Nhits2=Nhits2-1;
    }
    G4String proName = postStepPoint->GetProcessDefinedStep()->GetProcessName();
    // G4ProcessType proType = postStepPoint->GetProcessDefinedStep()->GetProcessType();
    // G4int proSubType = postStepPoint->GetProcessDefinedStep()->GetProcessSubType();
    
    //auto analysisManager = G4AnalysisManager::Instance();
    G4double engDep = aStep->GetTotalEnergyDeposit();
    if(engDep!=0){
        if(name=="opticalphoton"){
            if(PostVolume=="WorldVol_PV")Sum=Sum+1;
        }
    }
    if(name=="e-"){
        if(PreVolume=="CsIPhys"&&PostVolume=="ContainPhys")Eng_Ele_Lose=engDep+Eng_Ele_Lose;
        if(PreVolume=="ContainPhys"&&PostVolume=="ContainPhys")Eng_Ele_Lose=engDep+Eng_Ele_Lose;
        if(PreVolume=="ContainPhys"&&PostVolume=="BoxPhys")Eng_Ele_Lose=engDep+Eng_Ele_Lose;
        if(PreVolume=="BoxPhys"&&PostVolume=="BoxPhys")Eng_Ele_Lose=engDep+Eng_Ele_Lose;
        if(PreVolume=="BoxPhys"&&PostVolume=="WorldVol_PV")Eng_Ele_Lose=engDep+Eng_Ele_Lose;
    }
    if(parentID==0){
        if(engDep!=0)
        gamma_dep=engDep+gamma_dep;
    }
    if(name=="e-")e_eng=engDep+e_eng;

    if(parentID==0&&name=="gamma"){
        Eng_g_0_dop.push_back(aStep->GetTotalEnergyDeposit());
        Eng_g_0_dop_sum=engDep+Eng_g_0_dop_sum;
    }
    if(parentID==1&&name=="e-")Eng_e_1_dop.push_back(aStep->GetTotalEnergyDeposit());
    if(parentID==1&&name=="gamma")Eng_g_1_dop.push_back(aStep->GetTotalEnergyDeposit());
    //G4int pdgID = particle->GetPDGEncoding();
    //G4int parentID = aTrack->GetParentID();
   // if(engDep!=0 ){
        //Energy=engDep+Energy;
        //double Select=10;
        
     /*   if(name=="opticalphoton"){
                //SCintillation=SCintillation+1;
                Energy=engDep+Energy;
               // Eng_SC.push_back(aStep->GetTotalEnergyDeposit());
            }
        if(parentID==0){
            if(name=="e-")Select=0;//e-
            if(name=="gamma")Select=1;//
        }
        G4int idx = aTrack->GetCreatorModelID();
        if(idx==fPhotoAuger){
            Energy=engDep+Energy;
        }//Eng_Pho.push_back(aTrack->GetTotalEnergy());
        if(Select==1){
            if(name=="gamma"){
                
                //G4cout<<"gamma phy: "<<proName<<G4endl;
                //Energy=engDep+Energy;
            }
            if(name=="e-"){
                //G4cout<<"e- phy: "<<proName<<G4endl;
                if(proName=="eIoni"){
                    EIoni=EIoni+1;
                    Eng_Ei.push_back(aStep->GetTotalEnergyDeposit());
                }
                if(proName=="msc"){
                MSC=MSC+1;
                Eng_MSC.push_back(aStep->GetTotalEnergyDeposit());
                }
            }
            
            
        }
        if(Select==0){
            if(name=="e-"){
               // G4cout<<"phy: "<<proName<<G4endl;
                Energy=engDep+Energy;
                if(proName=="msc"){
                MSC=MSC+1;
                Eng_MSC.push_back(aStep->GetTotalEnergyDeposit());
                }
                else if(proName=="eIoni"){
                EIoni=EIoni+1;
                Eng_Ei.push_back(aStep->GetTotalEnergyDeposit());
                }
                else if(proName=="Scintillation"){
                SCintillation=SCintillation+1;
                Eng_SC.push_back(aStep->GetTotalEnergyDeposit());
                }
            }
            
        //if(name=="e-"){
            
            //if(parentID==0){
                E_num=E_num+1;
            //G4cout<<"proName:"<<proName<<G4endl;
        

        }*/
        //G4String creat_processName=aTrack->GetCreatorModelName();
        //G4cout<<"ModelName"<<creat_processName<<G4endl;
        
    //}
   // G4cout<<"the Vpre is"<< PostVolume <<G4endl;
    //G4cout<<"the Vpost is"<< PreVolume <<G4endl;
    //G4cout<<"the V2 is "<< thisVolume <<G4endl;
    //G4cout<<""<<G4endl;
    //G4cout<<"the V2 is"<< thisVolume2<<G4endl;
    //以下是G4Step常见的一些参数获取方法
    {
        /*
        G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
        G4StepPoint *postStepPoint = aStep->GetPostStepPoint();

        G4double stepLen = aStep->GetStepLength(); //the same as aTrack->GetTrackLength()
        G4double engDep = aStep->GetTotalEnergyDeposit(); // energy deposit inside this step
        G4double engDep2 = aStep->GetNonIonizingEnergyDeposit();
        G4bool isFirst = aStep->IsFirstStepInVolume();
        G4bool isLast = aStep->IsLastStepInVolume();
        G4cout << "-->" << stepLen << " " << engDep << " " << engDep2 << " " << isFirst << " " << isLast << G4endl;

        G4ThreeVector prePos = preStepPoint->GetPosition();
        G4double preGtime = preStepPoint->GetGlobalTime();
        G4double preLtime = preStepPoint->GetLocalTime();
        G4double prePtime = preStepPoint->GetProperTime();
        G4cout << "-->prePos=" << prePos << " " << preGtime << " " << preLtime << " " << prePtime << G4endl;
        G4ThreeVector postPos = postStepPoint->GetPosition(); // same as the information get from aTrack
        G4double postGtime = postStepPoint->GetGlobalTime();
        G4double postLtime = postStepPoint->GetLocalTime();
        G4double postPtime = postStepPoint->GetProperTime();
        G4cout << "-->postPos=" << postPos << " " << postGtime << " " << postLtime << " " << postPtime << G4endl;
        
        G4ThreeVector preMomDir = preStepPoint->GetMomentumDirection();
        G4ThreeVector preMomtum = preStepPoint->GetMomentum();
        G4double preEngTot = preStepPoint->GetTotalEnergy();
        G4double preEng = preStepPoint->GetKineticEnergy();
        G4double preV = preStepPoint->GetVelocity();
        G4double preBeta = preStepPoint->GetBeta();
        G4double preGamma = preStepPoint->GetGamma();

        G4ThreeVector postMomDir = postStepPoint->GetMomentumDirection();
        G4ThreeVector postMomtum = postStepPoint->GetMomentum();
        G4double postEngTot = postStepPoint->GetTotalEnergy();
        G4double postEng = postStepPoint->GetKineticEnergy();
        G4double postV = postStepPoint->GetVelocity();
        G4double postBeta = postStepPoint->GetBeta();
        G4double postGamma = postStepPoint->GetGamma();
        
        G4double mass = preStepPoint->GetMass();
        G4double charge = preStepPoint->GetCharge();
        G4double magmom = preStepPoint->GetMagneticMoment();

        auto *pVolume = postStepPoint->GetTouchableHandle()->GetVolume();
        auto *pLogicVolume = (pVolume!=NULL) ? pVolume->GetLogicalVolume() : NULL;

        */
    }

    //1.3 以下是拿DetectorConstruction相关参数的方法，根据具体情况修改即可
    {
        /*
        //拿Detector指针
        const MyDetectorConstruction* detectorConstruction = static_cast<const MyDetectorConstruction*> (G4RunManager::GetRunManager()->GetUserDetectorConstruction());

        //查找已知名字的逻辑体
        G4Box *fBox; 
        G4LogicalVolume* fLogicVol = G4LogicalVolumeStore::GetInstance()->GetVolume("YourLogicVolumeName");
        if ( fLogicVol ) fBox = dynamic_cast<G4Box*>(fLogicVol->GetSolid());
        fBox->GetXHalfLength() //这样就可以拿到具体的参数了。根据具体情况修改即可
        */
    }

    //-------
    //#ANALYSIS 4.4 在Steppinging的时候保存相应数据

    //---
    //2. 添加一些判断，并保存对应的数据。以下为演示，且按ANALYSIS 2. Ntuple定义的结构进行保存

    // G4ThreeVector postPos = postStepPoint->GetPosition();
    // G4ThreeVector postMomtum = postStepPoint->GetMomentum();

   
    //analysisManager->AddNtupleRow(1);
    //Ntuple2: 保存中性粒子信息：
    /*if (charge == 0 && (parentID == 1 && pdgID == 0)) //要求来自入射粒子，且是光子
    {
        if (proName != "Cerenkov") //只要切伦科夫过程
            return;

        auto *pVolume = postStepPoint->GetTouchableHandle()->GetVolume();
        if (pVolume == NULL)
            return;

        G4LogicalVolume *presentVolume = pVolume->GetLogicalVolume();
        if (presentVolume->GetName() != "FR4BoxVol") //只要光子打到阳极板上的情况
            return;

        G4double optEng = 0.0012398 / aTrack->GetKineticEnergy(); //convert to [nm]
        G4double optX = postPos.x();
        G4double optY = postPos.y();
        G4double optZ = postPos.z();

        // auto analysisManager = G4AnalysisManager::Instance();
        // analysisManager->FillNtupleDColumn(2, 0, optEng);
        // analysisManager->FillNtupleDColumn(2, 1, optX);
        // analysisManager->FillNtupleDColumn(2, 2, optY);
        // analysisManager->FillNtupleDColumn(2, 3, optZ);
        // analysisManager->FillNtupleDColumn(2, 4, aTrack->GetVertexPosition().x());
        // analysisManager->FillNtupleDColumn(2, 5, aTrack->GetVertexPosition().y());
        // analysisManager->FillNtupleDColumn(2, 6, aTrack->GetVertexPosition().z());
        // analysisManager->AddNtupleRow(2);
    }*/

    //Ntuple3: a)TrkLen在这里累加后，在EndOfEvent处保存； b)用vector来保存保存入射粒子的hit信息
 /*   if (parentID == 0)
    {
        G4double stepLen = aStep->GetStepLength();
        fTrkLen += stepLen;

        fEdeps.push_back(aStep->GetTotalEnergyDeposit());
        fHitsX.push_back(postPos.x());
        fHitsY.push_back(postPos.y());
        fHitsZ.push_back(postPos.z());
    }*/

    return;
}
